### Does Public Opinion on Foreign Policy Affect Elite Preferences?
### Evidence from the 2022 US Sanctions against Russia

### Anton Peez & Felix Bethke
### International Studies Quarterly

## @knitr Fig-2

rm(list = ls())

# Load/install all packages contained in the script
file_path <- "Public-TRIP-Code.R"
script <- readLines(file_path)

# Extract packages loaded with library() or require()
package_lines <- grep("^(library|require)\\(", script, value = TRUE)
package_names <- gsub(".*(library|require)\\(([^)]+)\\).*", "\\2", package_lines)
package_names <- gsub("[\"']", "", package_names)
package_names <- unique(package_names)

# Automatically install and load the packages using pacman
if (!requireNamespace("pacman", quietly = TRUE)) {
  install.packages("pacman")  # Install pacman if not already installed
}

# See the end of the script for exact information on package versions session details

### Script

library(tidyverse)
library(ggplot2)
library(ggthemes)
library(Rmisc)
library(here)
library(conflicted)
library(kableExtra)

conflict_prefer("rename", "dplyr")
conflict_prefer("here", "here")
conflict_prefer("filter", "dplyr")
conflict_prefer("summarize", "dplyr")

### Opening and summarizing the survey data

library(readr)
TRIP_Policymaker <- read_csv("TRIP_PolicymakerSurvey2022_Bethke_Peez[2].csv")

# Simple crosstab

Crosstab_Policymakers <- TRIP_Policymaker %>% group_by(prop_sanc) %>%
  dplyr::summarize(n = n(),
            increase         = sum(prop_sanc_q == "Increase"),
            unchanged        = sum(prop_sanc_q == "Continue unchanged"),
            decrease         = sum(prop_sanc_q == "Decrease"),
            share_increase   = increase/n  ,
            share_unchanged  = unchanged/n ,
            share_decrease   = decrease/n)

Crosstab_Policymakers <- TRIP_Policymaker %>% group_by(prop_sanc) %>%
  dplyr::mutate(n_status = n()) %>%
  group_by(prop_sanc_q, prop_sanc, n_status) %>%
  dplyr::summarize(n     = n()) 

Crosstab_Policymakers$share <- Crosstab_Policymakers$n/Crosstab_Policymakers$n_status













# Figure 1
# Confidence intervals

# Calculate the confidence intervals
Crosstab_Policymakers <- Crosstab_Policymakers %>%
  dplyr::mutate(ci_lower_95 = share - 1.96 * sqrt((share * (1 - share)) / n),
                ci_upper_95 = share + 1.96 * sqrt((share * (1 - share)) / n),
                ci_lower_90 = share - 1.645 * sqrt((share * (1 - share)) / n),
                ci_upper_90 = share + 1.645 * sqrt((share * (1 - share)) / n))

# Ensure the CIs don't go beyond the range of 0 to 1
Crosstab_Policymakers <- Crosstab_Policymakers %>%
  dplyr::mutate(ci_lower_95 = ifelse(ci_lower_95 < 0, 0, ci_lower_95),
                ci_upper_95 = ifelse(ci_upper_95 > 1, 1, ci_upper_95),
                ci_lower_90 = ifelse(ci_lower_90 < 0, 0, ci_lower_90),
                ci_upper_90 = ifelse(ci_upper_90 > 1, 1, ci_upper_90))

# Order the levels
Crosstab_Policymakers$prop_sanc_q <- factor(Crosstab_Policymakers$prop_sanc_q,
                                            levels = c("Increase",
                                                       "Continue unchanged",
                                                       "Decrease"))

Crosstab_Policymakers <- Crosstab_Policymakers %>% dplyr::arrange(prop_sanc_q)

# Rename to capitalize the labels
Crosstab_Policymakers$prop_sanc <- ifelse(Crosstab_Policymakers$prop_sanc == "control",
                                          "Control",
                                          Crosstab_Policymakers$prop_sanc)
Crosstab_Policymakers$prop_sanc <- ifelse(Crosstab_Policymakers$prop_sanc == "treatment",
                                          "Treatment",
                                          Crosstab_Policymakers$prop_sanc)


# Create a layer that will later be used for figure annotation
ann_text_0 <- data.frame(share   = 0.95,
                         prop_sanc = 1.5,
                         lab = "Text",
                         prop_sanc_q = factor("Increase", levels = c("Increase",
                                                                     "Continue unchanged",
                                                                     "Decrease")))

# Compile the figure
Fig_2 <- ggplot(Crosstab_Policymakers,
                aes(x = prop_sanc, y = share)) +
  geom_bar(stat = "identity", fill = "black", color = "black", width = 0.65) +
  facet_grid(. ~ prop_sanc_q) +
  theme_clean() +
  scale_fill_grey() +
  scale_y_continuous(limits = c(0.0, 1.0),
                     breaks = c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0),
                     labels = c("0%", " ", "20%", " ", "40%", " ", "60%", " ", "80%", " ", "100%")) +
  xlab(element_blank()) +
  ylab("Share supporting increasing sanctions") +
  labs(caption = "Two-sample test for equality of proportions (one-sided).\n90% and 95% confidence intervals indicated in grey.") +

  geom_errorbar(aes(ymin = ci_lower_95, ymax = ci_upper_95), 
                width = 0.125,  # Adjust the width of the whiskers
                color = "darkgrey", size = 0.8) +  # Thicker line for 95% CI whiskers
  
  geom_errorbar(aes(ymin = ci_lower_90, ymax = ci_upper_90), 
                width = 0.2,  # Adjust the width of the whiskers
                color = "darkgrey", size = 0.8) +   # Slightly thinner line for 90% CI whiskers
  
  geom_text(data = ann_text_0, label = "+8.3 pp., p=0.070")+
  theme(plot.background = element_rect(color = "black", fill = NA, size = 1.00))
  
Fig_2

ggsave("Fig_1.png", width = 6*1.3, height = 4.0)








### Main analysis -- Table 2
## @knitr Tab-1

# Subset -- increase/unchanged/decrease
# _i, _u, _d indicate increase, unchanged, decrease
xtab_i <- Crosstab_Policymakers %>%
  filter(prop_sanc_q == "Increase") %>%
  dplyr::arrange(desc(prop_sanc == "Treatment"))

xtab_u <- Crosstab_Policymakers %>%
  filter(prop_sanc_q == "Continue unchanged") %>%
  dplyr::arrange(desc(prop_sanc == "Treatment"))

xtab_d <- Crosstab_Policymakers %>%
  filter(prop_sanc_q == "Decrease") %>%
  dplyr::arrange(desc(prop_sanc == "Treatment"))

trip_prop_outcomes <- TRIP_Policymaker %>% group_by(prop_sanc, prop_sanc_q) %>% dplyr::summarise(n = n())
trip_prop_random   <- TRIP_Policymaker %>% group_by(prop_sanc) %>% dplyr::summarise(n = n())

# Increase -- one-/two-sided
inc_os        <- prop.test(x = xtab_i$n, n = xtab_i$n_status, alternative = "g", correct = F)
inc_os_90     <- prop.test(x = xtab_i$n, n = xtab_i$n_status, alternative = "g", correct = F, conf.level = 0.90)
inc_ts        <- prop.test(x = xtab_i$n, n = xtab_i$n_status, correct = F)

# Unchanged -- one-/two-sided
unc_os        <- prop.test(x = xtab_u$n, n = xtab_u$n_status, alternative = "g", correct = F)
unc_ts        <- prop.test(x = xtab_u$n, n = xtab_u$n_status, correct = F)

# Decrease -- one-/two-sided
dec_os <- prop.test(x = xtab_d$n, n = xtab_d$n_status, alternative = "g", correct = F)
dec_ts <- prop.test(x = xtab_d$n, n = xtab_d$n_status, correct = F)

# Format neatly
outcome    <- c("Increase", "Remain unchanged", "Decrease")
mean_treat    <- c(inc_os[["estimate"]][["prop 1"]],
                   unc_os[["estimate"]][["prop 1"]],
                   dec_os[["estimate"]][["prop 1"]])
mean_control  <- c(inc_os[["estimate"]][["prop 2"]],
                   unc_os[["estimate"]][["prop 2"]],
                   dec_os[["estimate"]][["prop 2"]])
p_value_os       <- c(as.numeric(inc_os[["p.value"]]),
                   "--",
                   "--") # No need for p-values as we have no expectations as to the direction
p_value_ts       <- c(as.numeric(inc_ts[["p.value"]]),
                      "--",
                      "--") # No need for p-values as we have no expectations as to the direction

# Filter for Treatment and Control groups
treatment <- Crosstab_Policymakers[Crosstab_Policymakers$prop_sanc == "Treatment", ]
control   <- Crosstab_Policymakers[Crosstab_Policymakers$prop_sanc == "Control", ]

# Create the words vectors
words_treat <- paste(treatment$n, "of", treatment$n_status)
words_control <- paste(control$n, "of", control$n_status)

# Bind
t_tests_main <- cbind(outcome, words_control, mean_control, words_treat, mean_treat, p_value_os, p_value_ts)
t_tests_main <- as.data.frame(t_tests_main)

t_tests_main$diff <-  as.numeric(t_tests_main$mean_treat) - as.numeric(t_tests_main$mean_control)

t_tests_main$mean_control <- round((as.numeric(t_tests_main$mean_control)  * 100), digits = 1)
t_tests_main$mean_treat   <- round((as.numeric(t_tests_main$mean_treat)    * 100), digits = 1)
t_tests_main$diff         <- round((as.numeric(t_tests_main$diff)          * 100), digits = 1)

t_tests_main$mean_control  <- paste0(as.character(t_tests_main$mean_control), "%")
t_tests_main$mean_treat    <- paste0(as.character(t_tests_main$mean_treat ), "%")
t_tests_main$diff          <- paste0(as.character(t_tests_main$diff ), " pp.")

t_tests_main$p_value_os      <- round(as.numeric(t_tests_main$p_value_os), digits = 3)
t_tests_main$p_value_os      <- ifelse(is.na(t_tests_main$p_value_os), "--", t_tests_main$p_value_os)

t_tests_main$p_value_ts      <- round(as.numeric(t_tests_main$p_value_ts), digits = 3)
t_tests_main$p_value_ts      <- ifelse(is.na(t_tests_main$p_value_ts), "--", t_tests_main$p_value_ts)

os_ci <- paste("[", round(inc_os[["conf.int"]][1]*100, digits = 1), "; ", "inf.", "]", sep = "")
ts_ci <- paste("[", round(inc_ts[["conf.int"]][1]*100, digits = 1), "; ", round(inc_ts[["conf.int"]][2]*100, digits = 1), "]", sep = "")

t_tests_main$ci95_os         <- c(os_ci, "--", "--")
t_tests_main$ci95_ts         <- c(ts_ci, "--", "--")

t_tests_main <- t_tests_main %>% relocate(diff,    .after = mean_treat)
t_tests_main <- t_tests_main %>% relocate(p_value_os, .after = ci95_os)
t_tests_main <- t_tests_main %>% relocate(p_value_ts, .after = ci95_ts)

t_tests_main_tab_os <- kbl(t_tests_main,
    caption = "Support for sanctions policies by treatment status.",
    booktabs = T,
    align = "lccccccccc",
    col.names = c("Policy preference",
                  "Respondents",
                  "Share",
                  "Respondents",
                  "Share",
                  "Diff.",
                  "95% CI",
                  "p",
                  "95% CI",
                  "p")) %>%
  kable_styling(font_size = 9,
                latex_options = "HOLD_position")         %>%
  add_header_above(c(" " = 1, "Control" = 2, "Treatment" = 2, "ATE" = 1, "one-sided" = 2, "two-sided" = 2)) %>%
  add_footnote("Significance tests are two-sample tests for equality of proportions.", notation="none")

t_tests_main_tab_os

readr::write_file(t_tests_main_tab_os, "Table_2.html")






## Retrospectively, given our effect size how many respondents would we have needed to reliably find it at p < 0.05, rather than p = 0.07?
#  (Footnote 18, Section 6.1)

p1    <- (83/122)     # Proportion in control group
p2    <- (100/131) # Proportion in treatment group (control + 8.3 percentage points)
alpha <- 0.05  # Significance level
power <- 0.80  # Desired power

# Perform power calculation
result <- power.prop.test(p1 = p1, p2 = p2, sig.level = alpha, power = power, alternative = "one.sided")

















### Table 1
#   Table comparing past similar experiments to ours

## @knitr Tab-2

# Read the menually-created spreadsheet summarizing other papers
library(readxl)
library(kableExtra)
options(scipen=999)

Studies <- read_excel("Studies.xlsx", col_types = c("text", 
                                                    "text", "text", "numeric", "text", "text", 
                                                    "text", "numeric"))

Studies$p <- ifelse(Studies$p < 0.001, "<0.001", Studies$p)
Studies$p <- ifelse(Studies$Study == "Present study",      "--",    Studies$p)

studies_tab <- kbl(Studies,
    caption = "Studies on elite responsiveness to public opinion on foreign policy.",
    booktabs = T,
    align = "lllcllcc",
    col.names = c("Study",
                  "Country",
                  "Elite sample",
                  "n",
                  "Policy decision",
                  "Hypotheticality",
                  "Effect",
                  "p")) %>%
  footnote(general = c("Significance tests are one-sided difference in proportions tests;",
           "calculated by this study's authors for Studies 2 and 3 for comparison across studies.")) %>%
  kable_styling(font_size = 7.5)

studies_tab

readr::write_file(studies_tab, "Table_1.html")


# The p-values in this table are derived as follows.

# Chu and Recchia
# Raw number are hand calculated based on appendix, but are slightly off from what their Table A2 shows
# They seem to do a two-sided diff in means test,
# But: (a) it should be one-sided based on theory (?), and (b) it should be difference in proportions (?)
prop_res_cease      <- prop.test(x = c(10,  9), n = c(52, 49), correct = F)
prop_res_continue_2 <- prop.test(x = c(33, 21), n = c(52, 49), alternative = "t", correct = F)
prop_res_continue_1 <- prop.test(x = c(33, 21), n = c(52, 49), alternative = "g", correct = F) # P-value shown in paper
prop_res_dko        <- prop.test(x = c(10, 18), n = c(52, 49), correct = F)

### Lin-Greenberg
#   Group sizes; Tab 5 in Lin-Greenberg's appendix: Public supports: N=100; public opposes: N=85
#   Support treatment: 97.6% favor strike; 0.976*100 = 97.60 favor strike
#   Oppose treatment:  83.0% favor strike; 0.830*85  = 70.55 favor strike
#   Likely some slight deviations bc of manipulation/attention checks, etc.
prop_strike_2 <- prop.test(x = c(97.60, 70.55), n = c(100, 85), alternative = "t", correct = F)
prop_strike_1 <- prop.test(x = c(97.60, 70.55), n = c(100, 85), alternative = "g", correct = F) # P-value shown in paper






# @knitr App-Tab-1

### Respondent descriptives

library(readr)
library(table1)
library(dplyr)
library(kableExtra)
library(forcats)

# Open the raw data again

demo <- read_csv("TRIP_PolicymakerSurvey2022_Bethke_Peez[2].csv")

# Rename and cluster some variables and categories

demo$Age <- fct_collapse(demo$Q276,
                         `25-44` = c("25-34", "35-44"))

demo <- demo %>%
  dplyr::mutate(Age = recode(Age,
                             "25-44" = '25--44',
                             "45-54" = '45--54',
                             "55-64" = '55--64',
                             "65-74" = "65--74",
                             "75-84" = "75--84" ))

demo$`Educational attainment` <- fct_collapse(demo$qg_108,
                                             `Professional (MPP/MPA/MIA, MBA, JD)` = c("JD", "MBA", "MPP/MPA/MIA"),
                                             `Other` = c("Other:", "All but dissertation"))

demo$Gender                   <- demo$qg_99
demo$`Economic views`         <- demo$qg_112b
demo$`Social views`           <- demo$qg_219b
demo$`Party identification`   <- demo$Q275
demo$`Job role`               <- demo$Q6
demo$Ethnicity                <- demo$qg_427

demo$Ethnicity <- ifelse((demo$Ethnicity == "White and non-white" | demo$Ethnicity == "Non-white"), "POC", demo$Ethnicity)

demo$`Years (list)`           <- demo$Year
demo$`Type of appointment`    <- demo$Q4

rank_desc <- demo %>%
  group_by(`Type of appointment` ) %>%
  dplyr::summarise(n = n())

demo$`Decisionmaking`    <- fct_collapse(demo$`Job role`,
                                       `Policy decision making` = c("Policy decision making"),
                                       `All other roles` = c("Support staff/other",
                                                             "Management",
                                                             "Policy analysis",
                                                             "Policy implementation",
                                                             "Policy support staff",
                                                             "Other:"))

demo$`Decisionmaking` <- factor(demo$`Decisionmaking`, levels = c("Policy decision making", "All other roles"))

decision_desc <- demo %>% group_by(`Decisionmaking`) %>% dplyr::summarise(n = n())

demo$`Job role`         <- fct_collapse(demo$Q6,
                                        `Support staff/other` = c("Other:", "Policy support staff"))

demo$`Expertise`        <- demo$Q95



# Re-order factors

demo <- demo %>% dplyr::mutate(`Educational attainment` = fct_relevel(`Educational attainment`, c("Bachelor's degree",
                                                           "Master's degree",
                                                           "Professional (MPP/MPA/MIA, MBA, JD)",
                                                           "Ph.D.",
                                                           "Other"
                                                           )))

demo <- demo %>% dplyr::mutate(`Ethnicity` = fct_relevel(`Ethnicity`, c(
                                                        "White",
                                                        "POC",
                                                        "No data"
                                                         )))

demo$`Party identification`         <- fct_collapse(demo$`Party identification`,
                                        `Independent/Other` = c("Independent", "Other:"))

demo <- demo %>% dplyr::mutate(`Party identification` = fct_relevel(`Party identification`, c(
                                                        "Republican",
                                                        "Independent/Other",
                                                        "Democrat"
                                                        )))

demo <- demo %>% dplyr::mutate(Expertise = fct_relevel(Expertise, c(
                                                                    "Security",
                                                                    "Development",
                                                                    "Trade"
                                                                  )))

demo$Gender <- ifelse(demo$Gender == "Prefer not to say", NA, demo$Gender)

# Rename control/treatment
demo$prop_sanc <- ifelse(demo$prop_sanc == "control",   "Control",   demo$prop_sanc)
demo$prop_sanc <- ifelse(demo$prop_sanc == "treatment", "Treatment", demo$prop_sanc)

levels(demo$`Party identification`) <- c("Republican",
                                         "Independent/Other",
                                         "Democrat")

# Function to calculate p-values for balance

pvalue <- function(x, ...) {
  # Construct vectors of data y, and groups (strata) g
  y <- unlist(x)
  g <- factor(rep(1:length(x), times=sapply(x, length)))
  if (is.numeric(y)) {
    # For numeric variables, perform a standard 2-sample t-test
    p <- t.test(y ~ g)$p.value
  } else {
    # For categorical variables, perform a chi-squared test of independence
    p <- chisq.test(table(y, g))$p.value
  }
  # Format the p-value, using an HTML entity for the less-than sign.
  # The initial empty string places the output on the line below the variable label.
  c("", sub("<", "&lt;", format.pval(p, digits=3, eps=0.001)))
}

## Balance

Tab_11a <- table1(~ 
                   Gender +
                   Age +
                   Ethnicity +
                    `Educational attainment` +
                    `Expertise`          +
                    `Decisionmaking`     +
                    `Party identification` | prop_sanc,    
                 data = demo,
                 overall=F,
                 extra.col=list(`p-value`=pvalue),
                 rowlabelhead = "Covariates",
                 caption = "Balance of covariates across treatment conditions.",
                 footnote = c("`Decisionmaking' clusters the `job task' categories as used in the subsequent analyses.", "`POC' clusters the following categories covered by TRIP: Asian Indian,", "Black/African-American, Hispanic and/or Latino, Japanese, Korean, Other."))

## Overall

Tab_11b <- table1(~ 
                    Gender +
                    Age +
                    Ethnicity +
                    `Educational attainment` +
                    `Expertise`          +
                    `Job role`     +
                    `Party identification`  ,    
                  data = demo,
                  rowlabelhead = "Covariates",
                  caption = "Descriptive statistics of the entire sample.",
                  footnote = c("`POC' clusters the following categories covered by TRIP: Asian Indian,", "Black/African-American, Hispanic and/or Latino, Japanese, Korean, Other."))

Table_A1 <- t1kable(Tab_11b,
        booktabs = T,
        longtable = T,
        align     = ('lrrr'),
        digits = c(NA, NA, NA, 3)
) %>%
  kable_styling(font_size = 10,
                latex_options = c("HOLD_position")) %>%
  kable_styling(latex_options = c("repeat_header"))

Table_A1

readr::write_file(Table_A1, "Table_A1.html")

# @knitr App-Tab-2

Table_A2 <- t1kable(Tab_11a,
        booktabs = T,
        longtable = T,
        align     = ('lrrr'),
        digits = c(NA, NA, NA, 3)
) %>%
  kable_styling(font_size = 10,
                latex_options = c("HOLD_position")) %>%
  kable_styling(latex_options = c("repeat_header"))

Table_A2

readr::write_file(Table_A2, "Table_A2.html")








# @knitr App-Tab-3
# Regression analyses -- Appendix

library(tidyverse)
library(modelsummary)
library(kableExtra)
library(estimatr)

# Some renaming
demo$partyid <- demo$`Party identification`
demo$edattn  <- demo$`Educational attainment`
  
# Binary DV
demo$prop_sanc_q_bin <- ifelse(demo$prop_sanc_q == "Increase", 1, 0)

logit_m1 <- glm(prop_sanc_q_bin ~ 
               prop_sanc +
               Age +
               Gender +
               Ethnicity +
               edattn +
               Expertise    +
               Decisionmaking     +
               partyid 
               ,
             family = "binomial", 
             demo)

#summary(logit_m1)

logit_m2 <- glm(prop_sanc_q_bin ~ 
                  prop_sanc 
                ,
                family = "binomial", 
                demo)

#summary(logit_m2)

# OLS
ols_m1 <- lm_robust(prop_sanc_q_bin ~ 
                      prop_sanc +
                      Age +
                      Gender +
                      Ethnicity +
                      edattn +
                      #ec_views +
                      #soc_views +
                      #`Years (list)` +
                      #`Type of appointment`  + 
                      #`Years in gov't (count)`+
                      Expertise +
                      Decisionmaking     +
                      partyid
              ,
              se_type = "stata",
              demo)

#summary(ols_m1)

ols_m2 <- lm_robust((prop_sanc_q_bin) ~ 
                      prop_sanc,
                se_type = "stata",
                data = demo)

#summary(ols_m2)

# Limit dataset to complete observations -- reviewer and editor comment

demo_lim <- demo %>% filter(  !is.na(Age)    & !is.na(edattn)
                            & !is.na(Gender)  # &!is.na(ec_views)
                            #& !is.na(soc_views) & !is.na(partyid)
                            & !is.na(Ethnicity) 
                            #& !is.na(`Years (list)`) 
                            #& !is.na(`Type of appointment`) 
                            #& !is.na(`Years in gov't (count)`) 
                            & !is.na(`Decisionmaking`) 
                            & !is.na(`Expertise`)
                            & !is.na(partyid))

ols_m3 <- lm_robust((prop_sanc_q_bin) ~ 
                      prop_sanc,
                   se_type = "stata",
                    data = demo_lim)

# summary(ols_m3)

# Principled regression using only relevant variables, as a reviewer sugegsted

ols_m4 <- lm_robust(prop_sanc_q_bin ~ 
                      prop_sanc +
                      #Age +
                      #edattn +
                      Gender +
                      #ec_views +
                      #soc_views +
                      #Ethnicity +
                      #`Years (list)` +
                      #`Type of appointment`  + 
                      #`Years in gov't (count)`+
                      `Expertise` +
                      `Decisionmaking`     +
                      partyid
                    ,
                    se_type = "stata",
                    demo)

#summary(ols_m4)

### Intermission: one-sided tests done above; but with this restricted sample.
desc_demo_lim <- demo_lim %>% group_by(prop_sanc, prop_sanc_q_bin) %>% summarize(n = n())
prop_res_i_1 <- prop.test(x = c(61, 83), n = c((33+61), (26+83)), alternative = "l", correct = F)
# The effect is indeed larger among respondents who provided demographic data.

### Format as a single regression table

models <- list(
  "M1a"                      = ols_m2,
  "M2"                       = ols_m1,
  "M1b"                      = ols_m3,
  "M3"                       = ols_m4
)

var_names <- c('(Intercept)'    = 'Intercept',
               'prop_sancTreatment' = 'Treated (ref.: control)',
               
               'GenderMale' = 'Gender identity: Male (ref.: female)',     ### Demographics
               'GenderPrefer not to say' = 'Gender identity: NA',
               
               'Age45--54' = 'Age: 45--54 (ref.: 25--44)',
               'Age55--64' = 'Age: 55--64',
               'Age65--74' = 'Age: 65--74',
               'Age75--84' = 'Age: 75--84',
               
               "EthnicityPOC" = "Ethnicity: POC (ref.: white)",
               "EthnicityNo data" = "Ethnicity: NA",
               
               "edattnMaster's degree" = "Education: Master's (ref.: Bachelor's)",
               "edattnProfessional (MPP/MPA/MIA, MBA, JD)" = "Education: Professional degree",
               "edattnPh.D." = "Education: Ph.D.",
               "edattnOther" = "Education: Other",
               
               "ExpertiseDevelopment" = "Expertise: Development (ref.: Security)",
               "ExpertiseTrade"    = "Expertise: Trade",
               
               "DecisionmakingPolicy decision making" = "Main job role: Decisionmaking",
               "DecisionmakingAll other roles" = "Main job role: All other tasks (ref.: Decisionmaking)",
               
               "partyidIndependent/Other" = "Party identification: Independent/Other (ref.: Republican)",   ### Political views
               "partyidDemocrat" = "Party identification: Democrat"
)

models_4 <- modelsummary::modelsummary(models,
                                       output = "kableExtra",
             longtable = T,                          
             coef_map = var_names,
             stars = c("†" = 0.1, "*" = 0.05, "**" = 0.01, "***" = 0.001),
             title = "OLS regression models with and without covariates.",
             gof_map = c("nobs", "r.squared"),
             #booktabs = T,
             notes = c("Models are OLS regressions with robust standard errors.", "``POC' clusters the following categories covered by TRIP: Asian Indian,", "Black/African-American, Hispanic and/or Latino, Japanese, Korean, white and non-white, Other.")) %>%
             #kbl(format = "latex") %>%
             kableExtra::kable_styling(font_size = 9, latex_options = c("HOLD_position", "scale_down"))  %>%
             kableExtra::add_header_above(c(" " = 1, "Outcome: increase sanctions Y/N" = 4))
  
models_4

readr::write_file(models_4, "Table_A3.html")






### Responses by field day.

# @knitr App-Fig-1

demo <- demo %>%
  dplyr::mutate(sub_date = lubridate::as_date(EndDate))

responses_by_day <- demo %>%
  group_by(sub_date) %>%
  summarize(N         = n()) %>%
  dplyr::mutate(Share        = N / sum(N),
         `Cum. share` = cumsum(Share))

Fig_Days <- ggplot(demo, aes(x = sub_date)) + 
  
  geom_vline(xintercept=as.Date("2022-12-15"), linetype = "dashed", color = "darkgrey")+
  annotate("text", x = as.Date("2022-12-17"), y = 25,
           label = "15 December 2022:\nNew US sanctions", hjust = 0,  vjust = 0.5,
           lineheight = 0.75, size = 3) +
  
  geom_vline(xintercept=as.Date("2022-12-07"), linetype = "dashed", color = "darkgrey")+
  annotate("text", x = as.Date("2022-12-09"), y = 35,
           label = "07 December 2022:\n9th EU sanctions package", hjust = 0,  vjust = 0.5,
           lineheight = 0.75, size = 3) +
  
  geom_vline(xintercept=as.Date("2022-12-03"), linetype = "dashed", color = "darkgrey")+
  annotate("text", x = as.Date("2022-12-09"), y = 45,
           label = "03 December 2022:\nG7 crude oil price cap set to $60", hjust = 0,  vjust = 0.5,
           lineheight = 0.75, size = 3) +  
  
  geom_vline(xintercept=as.Date("2022-11-30"), linetype = "dashed", color = "darkgrey")+
  annotate("text", x = as.Date("2022-11-28"), y = 35,
           label = "30 November 2022:\nEU proposal for a\nwar crimes tribunal", hjust = 1, vjust = 0.5,
           lineheight = 0.75, size = 3) +
  
  geom_histogram(binwidth = 1, fill = "darkgrey", color = "black") +
  labs(y = "Number of responses")+
  scale_x_date(date_labels = "%d %B",
               date_breaks = "1 week",
               limits = as.Date(c("2022-11-01", "2023-01-19")))+
  scale_y_continuous(limits = c(0, 50))+
  theme_clean()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.title.x = element_blank())

Fig_Days

ggsave("Fig_A5.png", width = 6, height = 3)






# @knitr App-Tab-4
# Naive monthly variables

library(lubridate)

demo$sub_month <- as.character(month(demo$EndDate, label = T))

demo$sub_month <- factor(demo$sub_month,
                         levels = c("Nov",
                                    "Dec",
                                    "Jan"))

ols_time_m1 <- lm_robust(prop_sanc_q_bin ~ 
                         prop_sanc +
                         sub_month,
                        se_type = "stata",
                        data = demo)

# summary(ols_time_m1)

# Binary pre-/post-event variables
# One day before

demo$pre_post_EU_trib <- ifelse(demo$EndDate > as.Date("2022-11-29"), 1, 0)
demo$pre_post_G7_cap  <- ifelse(demo$EndDate > as.Date("2022-12-02"), 1, 0)
demo$pre_post_EU_sanc <- ifelse(demo$EndDate > as.Date("2022-12-06"), 1, 0)

demo$pre_post_US_sanc <- ifelse(demo$EndDate > as.Date("2022-12-14"), 1, 0)

ols_time_m2 <- lm_robust(prop_sanc_q_bin ~ 
                           prop_sanc +
                           pre_post_EU_trib +
                           pre_post_G7_cap +
                           pre_post_EU_sanc,
                         se_type = "stata",
                         data = demo)

# summary(ols_time_m2)

ols_time_m3 <- lm_robust(prop_sanc_q_bin ~ 
                           prop_sanc +
                           pre_post_US_sanc,
                         se_type = "stata",
                         data = demo)

# summary(ols_time_m3)

# Duration of response
demo$sub_duration <- as.numeric((demo$EndDate - demo$StartDate) + 1)

ols_time_m4 <- lm_robust(prop_sanc_q_bin ~ 
                           prop_sanc +
                           log(sub_duration),
                         se_type = "stata",
                         data = demo)

# summary(ols_time_m4)

models <- list(
  "Months"   = ols_time_m1,
  "Events"   = ols_time_m2,
  "US sanctions"   = ols_time_m3,
  "Duration" = ols_time_m4
)

var_names <- c('(Intercept)'    = 'Intercept',
               'prop_sancTreatment' = 'Treated (ref.: control)',
               "sub_monthDec" = "Month: December (ref.: November)",
               "sub_monthJan" = "Month: January",
               "pre_post_EU_trib" = "Pre/post EU's tribunal proposal",
               "pre_post_G7_cap"  = "Pre/post G7 oil cap",
               "pre_post_EU_sanc" = "Pre/post 9th EU sanctions package",
               "pre_post_US_sanc" = "Pre/post new US sanctions",
               "log(sub_duration)" = "Duration of response in days (log)"
               )

models_time <- modelsummary::modelsummary(models,
                                          output = "kableExtra",
                                       coef_map = var_names,
                                       stars = c("†" = 0.1, "*" = 0.05, "**" = 0.01, "***" = 0.001),
                                       title = "The effects of survey submission date on support for sanctions.",
                                       gof_map = c("nobs", "r.squared"),
                                       #booktabs = T,
                                       notes = c("Models are OLS regressions with robust standard errors.")) %>%
  #kbl(format = "latex") %>%
  kableExtra::kable_styling(font_size = 9,
                            latex_options = "HOLD_position"
  )  %>%
  add_header_above(c(" " = 1, "Outcome: increase sanctions Y/N" = 4))

models_time

readr::write_file(models_time, "Table_A10.html")





# @knitr App-Tab-5
# Descriptive table by month

means_time_month <- demo %>% dplyr::summarise(mean_nov = mean(ifelse(sub_month == "Nov", prop_sanc_q_bin, NA), na.rm = T),
                                        mean_dec = mean(ifelse(sub_month == "Dec", prop_sanc_q_bin, NA), na.rm = T),
                                        mean_jan = mean(ifelse(sub_month == "Jan", prop_sanc_q_bin, NA), na.rm = T))

means_time_month$mean_nov       <- round((means_time_month$mean_nov       * 100), digits = 1)
means_time_month$mean_dec       <- round((means_time_month$mean_dec       * 100), digits = 1)
means_time_month$mean_jan       <- round((means_time_month$mean_jan       * 100), digits = 1)

means_time_month$mean_nov       <- paste0(as.character(means_time_month$mean_nov       ), "\\%")
means_time_month$mean_dec       <- paste0(as.character(means_time_month$mean_dec       ), "\\%")
means_time_month$mean_jan       <- paste0(as.character(means_time_month$mean_jan       ), "\\%")

means_time_month_tab <- kbl(means_time_month,
                            booktabs = T,
                            align     = ('ccc'),
                            escape = FALSE,
                            caption = "Share of respondents who support increasing sanctions, by month.",
                            linesep = "",
                            col.names = c("November 2022 (n = 70)",
                                          "December 2022 (n = 142)",
                                          "January 2023 (n = 41)")) %>%
  kable_styling(font_size = 10,
                latex_options = "HOLD_position")

means_time_month_tab

readr::write_file(means_time_month_tab, "Table_A11.html")







# @knitr App-Tab-6
# Descriptive table pre/post event

means_time_event <- demo %>% dplyr::summarise(
  mean_trib_pre  = mean(ifelse(pre_post_EU_trib == 0, prop_sanc_q_bin, NA), na.rm = T),
  mean_trib_post = mean(ifelse(pre_post_EU_trib == 1, prop_sanc_q_bin, NA), na.rm = T),
  mean_cap_pre   = mean(ifelse(pre_post_G7_cap  == 0, prop_sanc_q_bin, NA), na.rm = T),
  mean_cap_post  = mean(ifelse(pre_post_G7_cap  == 1, prop_sanc_q_bin, NA), na.rm = T),
  mean_sanc_pre  = mean(ifelse(pre_post_EU_sanc == 0, prop_sanc_q_bin, NA), na.rm = T),
  mean_sanc_post = mean(ifelse(pre_post_EU_sanc == 1, prop_sanc_q_bin, NA), na.rm = T),
  mean_us_sanc_pre  = mean(ifelse(pre_post_US_sanc == 0, prop_sanc_q_bin, NA), na.rm = T),
  mean_us_sanc_post = mean(ifelse(pre_post_US_sanc == 1, prop_sanc_q_bin, NA), na.rm = T))

means_time_event$mean_trib_pre  <- round((means_time_event$mean_trib_pre  * 100), digits = 1)
means_time_event$mean_cap_pre   <- round((means_time_event$mean_cap_pre   * 100), digits = 1)
means_time_event$mean_sanc_pre  <- round((means_time_event$mean_sanc_pre  * 100), digits = 1)
means_time_event$mean_us_sanc_pre  <- round((means_time_event$mean_us_sanc_pre  * 100), digits = 1)

means_time_event$mean_trib_post <- round((means_time_event$mean_trib_post * 100), digits = 1)
means_time_event$mean_cap_post  <- round((means_time_event$mean_cap_post  * 100), digits = 1)
means_time_event$mean_sanc_post <- round((means_time_event$mean_sanc_post * 100), digits = 1)
means_time_event$mean_us_sanc_post  <- round((means_time_event$mean_us_sanc_post  * 100), digits = 1)

means_time_event$mean_trib_pre  <- paste0(as.character(means_time_event$mean_trib_pre  ), "\\%")
means_time_event$mean_cap_pre   <- paste0(as.character(means_time_event$mean_cap_pre   ), "\\%")
means_time_event$mean_sanc_pre  <- paste0(as.character(means_time_event$mean_sanc_pre  ), "\\%")
means_time_event$mean_us_sanc_pre  <- paste0(as.character(means_time_event$mean_us_sanc_pre  ), "\\%")

means_time_event$mean_trib_post <- paste0(as.character(means_time_event$mean_trib_post ), "\\%")
means_time_event$mean_cap_post  <- paste0(as.character(means_time_event$mean_cap_post  ), "\\%")
means_time_event$mean_sanc_post <- paste0(as.character(means_time_event$mean_sanc_post ), "\\%")
means_time_event$mean_us_sanc_post  <- paste0(as.character(means_time_event$mean_us_sanc_post  ), "\\%")

means_time_event_tab <- kbl(means_time_event,
              booktabs = T,
              align     = ('cccccccc'),
              escape = FALSE,
              caption = "Share of respondents who support increasing sanctions, before and after important events.",
              linesep = "",
              col.names = c("Pre (58)",
                            "Post (195)",
                            "Pre (76)",
                            "Post (177)",
                            "Pre (117)",
                            "Post (136)",
                            "Pre (158)",
                            "Post (95)")) %>%
  add_header_above(c("Tribunal proposal" = 2, "G7 oil cap" = 2, "9th EU sanctions" = 2, "New US sanctions" = 2)) %>%
  kable_styling(font_size = 10,
                latex_options = "HOLD_position") %>%
  add_footnote("Values in parentheses are the observations in each group.", notation="none")

means_time_event_tab

readr::write_file(means_time_event_tab, "Table_A12.html")





### Heterogeneous treatment effects by expertise
### @knitr App-Tab-exp

### Treatment effects by subject-matter expertise -- two-sided t test

demo_sec   <- demo %>% filter(Expertise == "Security")
demo_dev   <- demo %>% filter(Expertise == "Development")
demo_trade <- demo %>% filter(Expertise == "Trade")

treat_sec   <- t.test(prop_sanc_q_bin ~ prop_sanc, data = demo_sec  )
treat_dev   <- t.test(prop_sanc_q_bin ~ prop_sanc, data = demo_dev  )
treat_trade <- t.test(prop_sanc_q_bin ~ prop_sanc, data = demo_trade)

expertise    <- c("Security", "Development", "Trade")
sample       <- c(nrow(demo_sec), nrow(demo_dev), nrow(demo_trade))
mean_treat   <- c(treat_sec[["estimate"]][["mean in group Treatment"]],
                  treat_dev[["estimate"]][["mean in group Treatment"]],
                  treat_trade[["estimate"]][["mean in group Treatment"]])
mean_control <- c(treat_sec[["estimate"]][["mean in group Control"]],
                    treat_dev[["estimate"]][["mean in group Control"]],
                    treat_trade[["estimate"]][["mean in group Control"]])
p_value      <- c(treat_sec[["p.value"]],
                  treat_dev[["p.value"]],
                  treat_trade[["p.value"]])

t_tests_expertise <- cbind(expertise, sample, mean_control, mean_treat, p_value)
t_tests_expertise <- as.data.frame(t_tests_expertise)

t_tests_expertise$diff <-  as.numeric(t_tests_expertise$mean_treat) - as.numeric(t_tests_expertise$mean_control)

t_tests_expertise$mean_control <- round((as.numeric(t_tests_expertise$mean_control)  * 100), digits = 1)
t_tests_expertise$mean_treat   <- round((as.numeric(t_tests_expertise$mean_treat)    * 100), digits = 1)
t_tests_expertise$diff   <- round((as.numeric(t_tests_expertise$diff)    * 100), digits = 1)

t_tests_expertise$mean_control  <- paste0(as.character(t_tests_expertise$mean_control), "\\%")
t_tests_expertise$mean_treat    <- paste0(as.character(t_tests_expertise$mean_treat ), "\\%")
t_tests_expertise$diff    <- paste0(as.character(t_tests_expertise$diff ), " pp.")

t_tests_expertise$p_value      <- round(as.numeric(t_tests_expertise$p_value), digits = 3)

t_tests_expertise <- t_tests_expertise %>% relocate(diff, .after = mean_treat)

t_tests_expertise_tab <- kbl(t_tests_expertise,
                            booktabs = T,
                            align     = ('lcccc'),
                            escape = FALSE,
                            caption = "Share of respondents who support increasing sanctions, by subject-matter expertise.",
                            linesep = "",
                            col.names = c("Expertise",
                                          "N",
                                          "Control",
                                          "Treated",
                                          "Difference",
                                          "p-value")) %>%
  add_header_above(c(" " = 2, "Share supporting" = 3, " " = 1)) %>%
  kable_styling(font_size = 10,
                latex_options = "HOLD_position") %>%
  add_footnote("Tests are Welch Two Sample t-tests (two-sided).", notation="none")

 #t_tests_expertise_tab

# Alternate procedure: two-sided diff in proportions

demo_sec   <- demo %>% filter(Expertise == "Security")
demo_dev   <- demo %>% filter(Expertise == "Development")
demo_trade <- demo %>% filter(Expertise == "Trade")

demo_sec_prop    <- demo_sec   %>%
  group_by(prop_sanc) %>%
  dplyr::summarise(n = n(), inc = sum(prop_sanc_q_bin)) %>%
  dplyr::arrange(desc(prop_sanc == "Treatment"))

demo_dev_prop    <- demo_dev   %>%
  group_by(prop_sanc) %>%
  dplyr::summarise(n = n(), inc = sum(prop_sanc_q_bin)) %>%
  dplyr::arrange(desc(prop_sanc == "Treatment"))

demo_trade_prop  <- demo_trade %>%
  group_by(prop_sanc) %>%
  dplyr::summarise(n = n(), inc = sum(prop_sanc_q_bin)) %>%
  dplyr::arrange(desc(prop_sanc == "Treatment"))

sec_os   <- prop.test(x = demo_sec_prop$inc,   n = demo_sec_prop$n,   alternative = "t", correct = F)
dev_os   <- prop.test(x = demo_dev_prop$inc,   n = demo_dev_prop$n,   alternative = "t", correct = F)
trade_os <- prop.test(x = demo_trade_prop$inc, n = demo_trade_prop$n, alternative = "t", correct = F)

expertise    <- c("Security", "Development", "Trade")
sample       <- c(nrow(demo_sec), nrow(demo_dev), nrow(demo_trade))
mean_treat   <- c(sec_os[["estimate"]][["prop 2"]],
                  dev_os[["estimate"]][["prop 2"]],
                  trade_os[["estimate"]][["prop 2"]])
mean_control <- c(sec_os[["estimate"]][["prop 1"]],
                  dev_os[["estimate"]][["prop 1"]],
                  trade_os[["estimate"]][["prop 1"]])
p_value      <- c(sec_os[["p.value"]],
                  dev_os[["p.value"]],
                  trade_os[["p.value"]])

t_tests_expertise <- cbind(expertise, sample, mean_control, mean_treat, p_value)
t_tests_expertise <- as.data.frame(t_tests_expertise)

t_tests_expertise_raw <- t_tests_expertise
write_csv(t_tests_expertise_raw, "t_tests_expertise_raw.csv")

t_tests_expertise$diff <-  as.numeric(t_tests_expertise$mean_treat) - as.numeric(t_tests_expertise$mean_control)

t_tests_expertise$mean_control <- round((as.numeric(t_tests_expertise$mean_control)  * 100), digits = 1)
t_tests_expertise$mean_treat   <- round((as.numeric(t_tests_expertise$mean_treat)    * 100), digits = 1)
t_tests_expertise$diff   <- round((as.numeric(t_tests_expertise$diff)    * 100), digits = 1)

t_tests_expertise$mean_control  <- paste0(as.character(t_tests_expertise$mean_control), "\\%")
t_tests_expertise$mean_treat    <- paste0(as.character(t_tests_expertise$mean_treat ), "\\%")
t_tests_expertise$diff    <- paste0(as.character(t_tests_expertise$diff ), "\\%")

t_tests_expertise$p_value      <- round(as.numeric(t_tests_expertise$p_value), digits = 3)

t_tests_expertise <- t_tests_expertise %>% relocate(diff, .after = mean_treat)

t_tests_expertise_tab_os <- kbl(t_tests_expertise,
                             booktabs = T,
                             align     = ('lccccc'),
                             escape = FALSE,
                             caption = "Share of respondents who support increasing sanctions, by subject-matter expertise.",
                             linesep = "",
                             col.names = c("Expertise",
                                           "N",
                                           "Control",
                                           "Treated",
                                           "Difference (ATE)",
                                           "p-value")) %>%
  add_header_above(c(" " = 2, "Share supporting" = 3, " " = 1)) %>%
  kable_styling(font_size = 10,
                latex_options = "HOLD_position") %>%
  add_footnote("Two-sample tests for equality of proportions (two-sided).", notation="none")

t_tests_expertise_tab_os

readr::write_file(t_tests_expertise_tab_os, "Table_A4.html")

# Alternate procedure: Interactions
# @knitr App-Fig-exp-b

lm_exp <- lm_robust(prop_sanc_q_bin ~ prop_sanc*Expertise, demo)

lm_exp_robust <- lm_robust(prop_sanc_q_bin ~ 
                           prop_sanc * Expertise,
                         se_type = "stata",
                         data = demo)

# summary(lm_exp)
# summary(lm_exp_robust)

var_names <- c('(Intercept)'          = 'Intercept',
               'prop_sancTreatment'   = 'Treated (ref.: control)',
               "ExpertiseDevelopment" = "Expertise: Development (ref.: Security)",
               "ExpertiseTrade"       = "Expertise: Trade",
               "prop_sancTreatment:ExpertiseDevelopment" = "Treatment × Expertise: Development",
               "prop_sancTreatment:ExpertiseTrade"       = "Treatment × Expertise: Trade")

models <- list("M1"   = lm_exp_robust)

models_exp <- modelsummary::modelsummary(models,
                                          coef_map = var_names,
                                          stars = c("†" = 0.1, "*" = 0.05, "**" = 0.01, "***" = 0.001),
                                          title = "The effects of treatment and subject-matter expertise on support for sanctions.",
                                          gof_map = c("nobs", "r.squared"),
                                          #booktabs = T,
                                          notes = c("Models are OLS regressions with robust standard errors.")) #%>%
  #::kable_styling(font_size = 9, latex_options = "HOLD_position"
  #)  %>%
  #add_header_above(c(" " = 1, "Outcome: increase sanctions Y/N" = 1))

### Marginal effects plot for expertise, based on this
library(margins)

##estimate linear model
model <- lm(prop_sanc_q_bin ~ prop_sanc * Expertise, data = demo)

##calculate marginal effects of treatment at different levels of domain
m <- margins(model, variables = c("prop_sanc"), at = list(Expertise = c("Security","Trade","Development")))

# plot(m)

##put all numerical results into object p
p <-summary(m)

#plot results from p with ggplot
margins_expertise <- ggplot(p, aes(x = Expertise, y = AME)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower,
                    ymax = upper,
                    width = 0.25)) +
  scale_y_continuous(breaks = c(-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6),
                     limits = c(-0.4, 0.6)) +
  geom_hline(yintercept = 0)+
  xlab("Subject-matter expertise")+
  ylab("Average Marginal Effect\nof the Treatment")+
  theme_clean()

margins_expertise

ggsave("Fig_A1.png", width = 4.5, height = 2)









### Gender

### @knitr App-Tab-gender

# Alternate procedure: two-sided diff in proportions

demo_male   <- demo %>% filter(Gender == "Male")
demo_female <- demo %>% filter(Gender == "Female")

demo_male_prop    <- demo_male   %>%
  group_by(prop_sanc) %>%
  dplyr::summarise(n = n(), inc = sum(prop_sanc_q_bin)) %>%
  dplyr::arrange(desc(prop_sanc == "Treatment"))

demo_female_prop  <- demo_female %>%
  group_by(prop_sanc) %>%
  dplyr::summarise(n = n(), inc = sum(prop_sanc_q_bin)) %>%
  dplyr::arrange(desc(prop_sanc == "Treatment"))

male_os   <- prop.test(x = demo_male_prop$inc,   n = demo_male_prop$n,   alternative = "t", correct = F)
female_os <- prop.test(x = demo_female_prop$inc, n = demo_female_prop$n, alternative = "t", correct = F)

gender    <- c("Male", "Female")
sample       <- c(nrow(demo_male), nrow(demo_female))
mean_treat   <- c(male_os[["estimate"]][["prop 2"]],
                  female_os[["estimate"]][["prop 2"]])
mean_control <- c(male_os[["estimate"]][["prop 1"]],
                  female_os[["estimate"]][["prop 1"]])
p_value      <- c(male_os[["p.value"]],
                  female_os[["p.value"]])

t_tests_gender <- cbind(gender, sample, mean_control, mean_treat, p_value)
t_tests_gender <- as.data.frame(t_tests_gender)

t_tests_gender_raw <- t_tests_gender
write_csv(t_tests_gender_raw, "t_tests_gender_raw.csv")

t_tests_gender$diff <-  as.numeric(t_tests_gender$mean_treat) - as.numeric(t_tests_gender$mean_control)

t_tests_gender$mean_control <- round((as.numeric(t_tests_gender$mean_control)  * 100), digits = 1)
t_tests_gender$mean_treat   <- round((as.numeric(t_tests_gender$mean_treat)    * 100), digits = 1)
t_tests_gender$diff   <- round((as.numeric(t_tests_gender$diff)    * 100), digits = 1)

t_tests_gender$mean_control  <- paste0(as.character(t_tests_gender$mean_control), "\\%")
t_tests_gender$mean_treat    <- paste0(as.character(t_tests_gender$mean_treat ), "\\%")
t_tests_gender$diff    <- paste0(as.character(t_tests_gender$diff ), "\\%")

t_tests_gender$p_value      <- round(as.numeric(t_tests_gender$p_value), digits = 3)

t_tests_gender <- t_tests_gender %>% relocate(diff, .after = mean_treat)

t_tests_gender_tab_os <- kbl(t_tests_gender,
                                booktabs = T,
                                align     = ('lccccc'),
                                escape = FALSE,
                                caption = "Share of respondents who support increasing sanctions, by respondent gender.",
                                linesep = "",
                                col.names = c("Gender",
                                              "N",
                                              "Control",
                                              "Treated",
                                              "Difference (ATE)",
                                              "p-value")) %>%
  add_header_above(c(" " = 2, "Share supporting" = 3, " " = 1)) %>%
  kable_styling(font_size = 10,
                latex_options = "HOLD_position") %>%
  add_footnote("Two-sample tests for equality of proportions (two-sided).", notation="none")

t_tests_gender_tab_os

readr::write_file(t_tests_gender_tab_os, "Table_A7.html")



# Check Gender x Expertise

# @knitr App-Tab-gender-desctab

demo_2 <- demo %>% filter(!is.na(Gender) & !is.na(Expertise))

# Expertise by Gender

bal_gender <- table1(~ 
                       Expertise | Gender,    
                     data = demo_2,
                     overall=F,
                     # extra.col=list(`p-value`=pvalue),
                     rowlabelhead = "Covariates",
                     caption = "Distribution of gender across expertise domains.")

distrib_exp_gender <- t1kable(bal_gender,
        booktabs = T,
        #longtable = T,
        align     = ('lrrr')
) %>%
  kable_styling(font_size = 10,
                latex_options = c("HOLD_position"))

distrib_exp_gender

readr::write_file(distrib_exp_gender, "Table_A8.html")



# Gender by expertise
# @knitr App-Tab-gender-desctab-2

bal_gen_exp <- table1(~ 
                       Gender | Expertise,    
                     data = demo_2,
                     overall=F,
                     # extra.col=list(`p-value`=pvalue),
                     rowlabelhead = "Covariates",
                     caption = "Distribution of expertise domains across gender.")

distrib_gender_exp <- t1kable(bal_gen_exp,
        booktabs = T,
        #longtable = T,
        align     = ('lrrr')
) %>%
  kable_styling(font_size = 10,
                latex_options = c("HOLD_position"))

distrib_gender_exp

readr::write_file(distrib_gender_exp, "Table_A9.html")



# Alternate procedure: Interactions
# @knitr App-Fig-gender

lm_gender <- lm_robust(prop_sanc_q_bin ~ prop_sanc*Gender, demo)

lm_gender_robust <- lm_robust(prop_sanc_q_bin ~ 
                             prop_sanc * Gender,
                           se_type = "stata",
                           data = demo)

# summary(lm_exp)
# summary(lm_exp_robust)

var_names <- c('(Intercept)'          = 'Intercept',
               'prop_sancTreatment'   = 'Treated (ref.: control)',
               "GenderMale" = "Gender: Male (ref.: Female)",
               "prop_sancTreatment:GenderMale" = "Treatment × Gender: Male")

models <- list("M1"   = lm_gender_robust)

models_gender <- modelsummary::modelsummary(models,
                                         coef_map = var_names,
                                         stars = c("†" = 0.1, "*" = 0.05, "**" = 0.01, "***" = 0.001),
                                         title = "The effects of treatment and respondent gender on support for sanctions.",
                                         gof_map = c("nobs", "r.squared"),
                                         #booktabs = T,
                                         notes = c("Models are OLS regressions with robust standard errors.")) #%>%
#::kable_styling(font_size = 9, latex_options = "HOLD_position"
#)  %>%
#add_header_above(c(" " = 1, "Outcome: increase sanctions Y/N" = 1))


### Marginal effects plot for expertise, based on this
library(margins)

##estimate linear model
model <- lm(prop_sanc_q_bin ~ prop_sanc * Gender, data = demo)

##calculate marginal effects of treatment at different levels of domain
m <- margins(model, variables = c("prop_sanc"), at = list(Gender = c("Male", "Female")))

# plot(m)

##put all numerical results into object p
p <-summary(m)

#plot results from p with ggplot
margins_gender <- ggplot(p, aes(x = Gender, y = AME)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower,
                    ymax = upper,
                    width = 0.15)) +
  scale_y_continuous(breaks = c(-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6),
                     limits = c(-0.4, 0.6)) +
  geom_hline(yintercept = 0)+
  xlab("Gender")+
  ylab("Average Marginal Effect\nof the Treatment")+
  theme_clean()

margins_gender

ggsave("Fig_A4.png", width = 4.5, height = 2)







### Decisionmaking vs. all other roles -- "skin in the game"

#   @knitr App-Tab-dec

demo_dec           <- demo %>%
  filter(`Decisionmaking` == "Policy decision making")

demo_other_roles   <- demo %>%
  filter(`Decisionmaking` == "All other roles")

demo_dec_prop         <- demo_dec          %>%
  group_by(prop_sanc) %>%
  dplyr::summarise(n = n(), inc = sum(prop_sanc_q_bin)) %>%
  dplyr::arrange(desc(prop_sanc == "Treatment"))

demo_other_roles_prop <- demo_other_roles  %>%
  group_by(prop_sanc) %>%
  dplyr::summarise(n = n(), inc = sum(prop_sanc_q_bin)) %>%
  dplyr::arrange(desc(prop_sanc == "Treatment"))

dec_os           <- prop.test(x = demo_dec_prop$inc,         n = demo_dec_prop$n,         alternative = "t", correct = F)
other_roles_os   <- prop.test(x = demo_other_roles_prop$inc, n = demo_other_roles_prop$n, alternative = "t", correct = F)

decision    <- c("Policy decisionmaking", "All other roles")
sample       <- c(nrow(demo_dec), nrow(demo_other_roles))
mean_treat   <- c(dec_os[["estimate"]][["prop 2"]],
                  other_roles_os[["estimate"]][["prop 2"]])
mean_control <- c(dec_os[["estimate"]][["prop 1"]],
                  other_roles_os[["estimate"]][["prop 1"]])
p_value      <- c(dec_os[["p.value"]],
                  other_roles_os[["p.value"]])

t_tests_dec <- cbind(decision, sample, mean_control, mean_treat, p_value)
t_tests_dec <- as.data.frame(t_tests_dec)

t_tests_dec_raw <- t_tests_dec
write_csv(t_tests_dec_raw, "t_tests_dec_raw.csv")

t_tests_dec$diff <-  as.numeric(t_tests_dec$mean_treat) - as.numeric(t_tests_dec$mean_control)

t_tests_dec$mean_control <- round((as.numeric(t_tests_dec$mean_control)  * 100), digits = 1)
t_tests_dec$mean_treat   <- round((as.numeric(t_tests_dec$mean_treat)    * 100), digits = 1)
t_tests_dec$diff   <- round((as.numeric(t_tests_dec$diff)    * 100), digits = 1)

t_tests_dec$mean_control  <- paste0(as.character(t_tests_dec$mean_control), "\\%")
t_tests_dec$mean_treat    <- paste0(as.character(t_tests_dec$mean_treat ), "\\%")
t_tests_dec$diff    <- paste0(as.character(t_tests_dec$diff ), " pp.")

t_tests_dec$p_value      <- round(as.numeric(t_tests_dec$p_value), digits = 3)

t_tests_dec <- t_tests_dec %>% relocate(diff, .after = mean_treat)

t_tests_dec_tab_os <- kbl(t_tests_dec,
                           booktabs = T,
                           align     = ('lccccc'),
                           escape = FALSE,
                           caption = "Share of respondents who support increasing sanctions, by main job task.",
                           linesep = "",
                           col.names = c("Main job task",
                                         "N",
                                         "Control",
                                         "Treated",
                                         "Difference (ATE)",
                                         "p-value")) %>%
  add_header_above(c(" " = 2, "Share supporting" = 3, " " = 1)) %>%
  kable_styling(font_size = 10,
                latex_options = "HOLD_position") %>%
  add_footnote("Two-sample tests for equality of proportions (two-sided).", notation="none")

t_tests_dec_tab_os

readr::write_file(t_tests_dec_tab_os, "Table_A5.html")

# Alternate procedure: Interactions
# @knitr App-Fig-dec-b

lm_dec <- lm_robust(prop_sanc_q_bin ~ prop_sanc*`Decisionmaking`, demo)

lm_dec_robust <- lm_robust(prop_sanc_q_bin ~ 
                              prop_sanc * `Decisionmaking`,
                            se_type = "stata",
                            data = demo)

# summary(lm_dec)
# summary(lm_dec_robust)

var_names <- c('(Intercept)'          = 'Intercept',
               'prop_sancTreatment'   = 'Treated (ref.: control)',
               "DecisionmakingPolicy decision making " = "Job role: Policy decisionmaking (ref.: all other roles)",
               "prop_sancTreatment:DecisionmakingPolicy decision making"         = "Treatment × Job role: Policy decisionmaking")

models <- list("M1"   = lm_dec_robust)

models_dec <- modelsummary::modelsummary(models,
                                         coef_map = var_names,
                                         stars = c("†" = 0.1, "*" = 0.05, "**" = 0.01, "***" = 0.001),
                                         title = "The effects of treatment and job responsbility/role on support for sanctions.",
                                         gof_map = c("nobs", "r.squared"),
                                         #booktabs = T,
                                         notes = c("Models are OLS regressions with robust standard errors.")) #%>%
  #::kable_styling(font_size = 9,
  #                          latex_options = "HOLD_position"
  #)  %>%
  #add_header_above(c(" " = 1, "Outcome: increase sanctions Y/N" = 1))

# models_dec

### Marginal effects plot for expertise, based on this
library(margins)

##estimate linear model
model <- lm(prop_sanc_q_bin ~ prop_sanc * `Decisionmaking`, data = demo)

##calculate marginal effects of treatment at different levels of domain
m <- margins(model, variables = c("prop_sanc"), at = list(`Decisionmaking` = c("Policy decision making", "All other roles")))

# plot(m)

##put all numerical results into object p
p <-summary(m)

#plot results from p with ggplot
margins_dec <- ggplot(p, aes(x = `Decisionmaking`, y = AME)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower,
                    ymax = upper,
                    width = 0.15)) +
  scale_y_continuous(breaks = c(-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6),
                     limits = c(-0.4, 0.6)) +
  geom_hline(yintercept = 0)+
  xlab("Primary job task")+
  ylab("Average Marginal Effect\nof the Treatment")+
  theme_clean()

margins_dec

ggsave("Fig_A2.png", width = 4.5, height = 2)






### Party identification
#   @knitr App-Tab-partyid

demo_dem           <- demo %>% filter(partyid == "Democrat")
demo_rep           <- demo %>% filter(partyid == "Republican")
demo_ind           <- demo %>% filter(partyid == "Independent/Other")

demo_dem_prop   <- demo_dem     %>%
  group_by(prop_sanc) %>%
  dplyr::summarise(n = n(), inc = sum(prop_sanc_q_bin)) 

demo_rep_prop   <- demo_rep     %>%
  group_by(prop_sanc) %>%
  dplyr::summarise(n = n(), inc = sum(prop_sanc_q_bin)) 

demo_ind_prop   <- demo_ind     %>%
  group_by(prop_sanc) %>%
  dplyr::summarise(n = n(), inc = sum(prop_sanc_q_bin)) 

dem_os   <- prop.test(x = demo_dem_prop$inc, n = demo_dem_prop$n, alternative = "t", correct = F)
rep_os   <- prop.test(x = demo_rep_prop$inc, n = demo_rep_prop$n, alternative = "t", correct = F)
ind_os   <- prop.test(x = demo_ind_prop$inc, n = demo_ind_prop$n, alternative = "t", correct = F)

partyid    <- c("Democrat", "Republican", "Independent/Other")
sample       <- c(nrow(demo_dem), nrow(demo_rep), nrow(demo_ind))
mean_treat   <- c(dem_os[["estimate"]][["prop 2"]],
                  rep_os[["estimate"]][["prop 2"]],
                  ind_os[["estimate"]][["prop 2"]])
mean_control <- c(dem_os[["estimate"]][["prop 1"]],
                  rep_os[["estimate"]][["prop 1"]],
                  ind_os[["estimate"]][["prop 1"]])
p_value      <- c(dem_os[["p.value"]],
                  rep_os[["p.value"]],
                  ind_os[["p.value"]])

t_tests_dem <- cbind(partyid, sample, mean_control, mean_treat, p_value)
t_tests_dem <- as.data.frame(t_tests_dem)

t_tests_dem_raw <- t_tests_dem
write_csv(t_tests_dem_raw, "t_tests_dem_raw.csv")

t_tests_dem$diff <-  as.numeric(t_tests_dem$mean_treat) - as.numeric(t_tests_dem$mean_control)

t_tests_dem$mean_control <- round((as.numeric(t_tests_dem$mean_control)  * 100), digits = 1)
t_tests_dem$mean_treat   <- round((as.numeric(t_tests_dem$mean_treat)    * 100), digits = 1)
t_tests_dem$diff   <- round((as.numeric(t_tests_dem$diff)    * 100), digits = 1)

t_tests_dem$mean_control  <- paste0(as.character(t_tests_dem$mean_control), "\\%")
t_tests_dem$mean_treat    <- paste0(as.character(t_tests_dem$mean_treat ), "\\%")
t_tests_dem$diff    <- paste0(as.character(t_tests_dem$diff ), " pp.")

t_tests_dem$p_value      <- round(as.numeric(t_tests_dem$p_value), digits = 3)

t_tests_dem <- t_tests_dem %>% relocate(diff, .after = mean_treat)

t_tests_dem_tab_os <- kbl(t_tests_dem,
                          booktabs = T,
                          align     = ('lccccc'),
                          escape = FALSE,
                          caption = "Share of respondents who support increasing sanctions, by party identification.",
                          linesep = "",
                          col.names = c("Party identification",
                                        "N",
                                        "Control",
                                        "Treated",
                                        "Difference (ATE)",
                                        "p-value")) %>%
  add_header_above(c(" " = 2, "Share supporting" = 3, " " = 1)) %>%
  kable_styling(font_size = 10,
                latex_options = "HOLD_position") %>%
  add_footnote("Two-sample tests for equality of proportions (two-sided).", notation="none")

t_tests_dem_tab_os

readr::write_file(t_tests_dem_tab_os, "Table_A6.html")

# Alternate procedure: Interactions
# @knitr App-Fig-partyid

lm_partyid <- lm_robust(prop_sanc_q_bin ~ prop_sanc*partyid, demo)

lm_partyid_robust <- lm_robust(prop_sanc_q_bin ~ 
                             prop_sanc * partyid,
                           se_type = "stata",
                           data = demo)

# summary(lm_partyid)
# summary(lm_partyid_robust)

### Marginal effects plot for expertise, based on this
library(margins)

##estimate linear model
model <- lm(prop_sanc_q_bin ~ prop_sanc * partyid, data = demo)

##calculate marginal effects of treatment at different levels of domain
m <- margins(model, variables = c("prop_sanc"), at = list(partyid = c("Democrat", "Republican", "Independent/Other")))

# plot(m)

##put all numerical results into object p
p <-summary(m)

#plot results from p with ggplot
margins_partyid <- ggplot(p, aes(x = partyid, y = AME)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower,
                    ymax = upper,
                    width = 0.15)) +
  scale_y_continuous(breaks = c(-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5),
                     limits = c(-0.4, 0.5)) +
  geom_hline(yintercept = 0)+
  xlab("Party identification")+
  ylab("Average Marginal Effect\nof the Treatment")+
  theme_clean()

margins_partyid

ggsave("Fig_A3.png", width = 4.5, height = 2)








### Figure in main text of all subgroup analyses
### @knitr fig-hetero

t_tests_expertise_raw <- read_csv("t_tests_expertise_raw.csv")
t_tests_gender_raw    <- read_csv("t_tests_gender_raw.csv")
t_tests_dec_raw       <- read_csv("t_tests_dec_raw.csv")
t_tests_dem_raw       <- read_csv("t_tests_dem_raw.csv")

t_tests_expertise_raw$subgroup <- "Expertise"
t_tests_gender_raw$subgroup    <- "Gender"
t_tests_dec_raw$subgroup       <- "Primary job task"
t_tests_dem_raw$subgroup       <- "Party identification"

t_tests_expertise_raw$ate <- t_tests_expertise_raw$mean_treat - t_tests_expertise_raw$mean_control
t_tests_gender_raw$ate      <- t_tests_gender_raw$mean_treat - t_tests_gender_raw$mean_control
t_tests_dec_raw$ate       <- t_tests_dec_raw$mean_treat - t_tests_dec_raw$mean_control
t_tests_dem_raw$ate       <- t_tests_dem_raw$mean_treat - t_tests_dem_raw$mean_control

wide_expertise <- t_tests_expertise_raw %>%
  pivot_longer(cols = starts_with("mean"),
               names_to = "status",
               values_to = "mean") %>%
  dplyr::mutate(status = if_else(status == "mean_treat", "Treatment", "Control")) %>%
  #dplyr::mutate(status = fct_relevel(status, c("Treatment","Control"))) %>%
  select(expertise, status, mean, everything())

wide_gender <- t_tests_gender_raw %>%
  pivot_longer(cols = starts_with("mean"),
               names_to = "status",
               values_to = "mean") %>%
  dplyr::mutate(status = if_else(status == "mean_treat", "Treatment", "Control")) %>%
  #dplyr::mutate(status = fct_relevel(status, c("Treatment","Control"))) %>%
  select(gender, status, mean, everything())

wide_task <- t_tests_dec_raw %>%
  pivot_longer(cols = starts_with("mean"),
               names_to = "status",
               values_to = "mean") %>%
  dplyr::mutate(status = if_else(status == "mean_treat", "Treatment", "Control")) %>%
  #dplyr::mutate(status = fct_relevel(status, c("Treatment","Control"))) %>%
  select(decision, status, mean, everything())

wide_partyid <- t_tests_dem_raw %>%
  pivot_longer(cols = starts_with("mean"),
               names_to = "status",
               values_to = "mean") %>%
  dplyr::mutate(status = if_else(status == "mean_treat", "Treatment", "Control")) %>%
  #dplyr::mutate(status = fct_relevel(status, c("Treatment","Control"))) %>%
  select(partyid, status, mean, everything())

# Add empty rows where needed

wide_task           <- rbind(wide_task,
                                  c(" ", "Control", NA, NA, NA, NA, NA))
wide_task           <- rbind(wide_task,
                                  c(" ", "Treatment", NA, NA, NA, NA, NA))

wide_gender           <- rbind(wide_gender,
                             c(" ", "Control", NA, NA, NA, NA, NA))
wide_gender           <- rbind(wide_gender,
                             c(" ", "Treatment", NA, NA, NA, NA, NA))

# Reorder the levels for each characteristic
new_order_expertise <- c("Development", "Security", "Trade", " ") 
new_order_gender      <- c("Male", "Female", " ")
new_order_task      <- c("Policy decisionmaking", "All other roles", " ", "  ")
new_order_partyid   <- c("Democrat", "Independent/Other", "Republican", " ")

wide_expertise$expertise <- factor(wide_expertise$expertise, levels = new_order_expertise)
wide_gender$gender       <- factor(wide_gender$gender, levels = new_order_gender)
wide_task$decision       <- factor(wide_task$decision, levels = new_order_task)
wide_partyid$partyid     <- factor(wide_partyid$partyid, levels = new_order_partyid)

## Custom annotations

ann_text_1 <- data.frame(mean   = 0.95,
                       status = 1.5,
                       lab = "Text",
                       expertise = factor("Development", levels = c("Development",
                                                                    "Security",
                                                                    "Trade",
                                                                    " ")))

ann_text_3 <- data.frame(mean   = 0.95,
                         status = 1.5,
                         lab = "Text",
                         decision = factor("Policy decisionmaking", levels = c("Policy decisionmaking",
                                                                               "All other roles")))

ann_text_4 <- data.frame(mean   = 0.95,
                         status = 1.5,
                         lab = "Text",
                         expertise = factor("Democrat", levels = c("Democrat",
                                                                   "Independent/Other",
                                                                   "Republican")))

ann_text_5 <- data.frame(mean   = 0.95,
                         status = 1.5,
                         lab = "Text",
                         gender = factor("Female",
                                       levels = c("Male",
                                                  "Female")))





## Figures

fig_exp <- ggplot(wide_expertise,
                        aes(x = status, y = as.numeric(mean))) +
  geom_bar(stat = "identity", fill = "black", color = "black", width = 0.65) +
  facet_grid(.~expertise) +
  theme_clean() +
  scale_fill_grey() +
  scale_y_continuous(limits = c(0.0, 1.0),
                     breaks = c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0),
                     labels = c("0%", " ", "20%", " ", "40%", " ", "60%", " ", "80%", " ", "100%")) +
  xlab(element_blank()) +
  ylab("Share supporting increasing sanctions") +
  geom_text(data = ann_text_1, label = "+24 pp., p=0.069") # Check!

gender_names <- c(
  `Male` = "Male",
  `Female` = "Female"
  )

fig_gender <- ggplot(wide_gender,
                  aes(x = status, y = as.numeric(mean))) +
  geom_bar(stat = "identity", fill = "black", color = "black", width = 0.65) +
  facet_grid(.~gender) +
  theme_clean() +
  scale_fill_grey() +
  scale_y_continuous(limits = c(0.0, 1.0),
                     breaks = c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0),
                     labels = c("0%", " ", "20%", " ", "40%", " ", "60%", " ", "80%", " ", "100%")) +
  xlab(element_blank()) +
  ylab("Share supporting increasing sanctions") +
  geom_text(data = ann_text_5, label = "+32.6 pp., p=0.023") +
  labs(caption = "Two-sample tests for equality of proportions (two-sided).\nResults with p < 0.10 shown.")

fig_task <- ggplot(wide_task,
                   aes(x = status, y = as.numeric(mean))) +
  geom_bar(stat = "identity", fill = "black", color = "black", width = 0.65) +
  facet_grid(.~decision) +
  theme_clean() +
  scale_fill_grey() +
  scale_y_continuous(limits = c(0.0, 1.0),
                     breaks = c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0),
                     labels = c("0%", " ", "20%", " ", "40%", " ", "60%", " ", "80%", " ", "100%")) +
  xlab(element_blank()) +
  ylab("Share supporting increasing sanctions") +
  geom_text(data = ann_text_3, label = "+17.9 pp., p=0.037")

fig_partyid <- ggplot(wide_partyid,
                   aes(x = status, y = as.numeric(mean))) +
  geom_bar(stat = "identity", fill = "black", color = "black", width = 0.65) +
  facet_grid(.~partyid) +
  theme_clean() +
  scale_fill_grey() +
  scale_y_continuous(limits = c(0.0, 1.0),
                     breaks = c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0),
                     labels = c("0%", " ", "20%", " ", "40%", " ", "60%", " ", "80%", " ", "100%")) +
  xlab(element_blank()) +
  ylab("Share supporting increasing sanctions") +
  geom_text(data = ann_text_4, label = c(" "))

## Combine

library(patchwork)

comb_hetero <- (fig_exp/fig_task/fig_partyid/fig_gender) + plot_annotation(tag_levels ='A')

comb_hetero

ggsave("Fig_2.png", width = 7.8, height = 11.424)






### Final step: Create a file with all session information
### The last few lines in this file indicate the versions used when writing and reproducing this script.

# Create a file to store session information
output_file <- "session_info.txt"

# Open a connection to the file
fileConn <- file(output_file)

# Capture session information with all loaded packages and versions
session_info <- capture.output(sessionInfo())

# Check if running in RStudio and get RStudio version
if (requireNamespace("rstudioapi", quietly = TRUE)) {
  rstudio_version <- rstudioapi::versionInfo()$version
  session_info <- c(session_info, "", paste("RStudio Version:", rstudio_version))
} else {
  session_info <- c(session_info, "", "RStudio Version: Not Available (not running in RStudio)")
}

# Get Base R version information
r_version <- paste("R Version:", R.version.string)
session_info <- c(session_info, r_version)

# Write all session info to the file
writeLines(session_info, fileConn)

# Close the file connection
close(fileConn)

# The exact session and package details for this code are as follows.

#R version 4.3.1 (2023-06-16)
#Platform: x86_64-apple-darwin20 (64-bit)
#Running under: macOS Sonoma 14.6.1
#
#Matrix products: default
#BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
#LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
#
#locale:
#  [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#
#time zone: Europe/Berlin
#tzcode source: internal
#
#attached base packages:
#  [1] stats     graphics  grDevices utils     datasets  methods   base     
#
#other attached packages:
#  [1] patchwork_1.2.0    margins_0.3.27     estimatr_1.0.4     modelsummary_2.1.1 table1_1.4.3       readxl_1.4.3       kableExtra_1.4.0
#[8] conflicted_1.2.0   here_1.0.1         Rmisc_1.5.1        plyr_1.8.9         lattice_0.22-6     ggthemes_5.1.0     lubridate_1.9.3   
#[15] forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4        purrr_1.0.2        readr_2.1.5        tidyr_1.3.1        tibble_3.2.1      
#[22] ggplot2_3.5.1      tidyverse_2.0.0   
#
#loaded via a namespace (and not attached):
#  [1] tidyselect_1.2.1    viridisLite_0.4.2   farver_2.1.2        fastmap_1.2.0       prediction_0.3.18   bayestestR_0.13.2  
#[7] pacman_0.5.1        digest_0.6.36       timechange_0.3.0    estimability_1.5.1  lifecycle_1.0.4     magrittr_2.0.3     
#[13] compiler_4.3.1      rlang_1.1.4         sass_0.4.9          tools_4.3.1         utf8_1.2.4          yaml_2.3.9         
#[19] data.table_1.15.4   knitr_1.48          labeling_0.4.3      bit_4.0.5           xml2_1.3.6          tinytable_0.3.0    
#[25] withr_3.0.0         datawizard_0.10.0   grid_4.3.1          fansi_1.0.6         xtable_1.8-4        colorspace_2.1-0   
#[31] future_1.33.2       MASS_7.3-60         emmeans_1.10.2      globals_0.16.3      scales_1.3.0        insight_0.20.1     
#[37] cli_3.6.3           mvtnorm_1.2-5       rmarkdown_2.27      crayon_1.5.3        ragg_1.3.2          generics_0.1.3     
#[43] performance_0.11.0  rstudioapi_0.16.0   future.apply_1.11.2 tzdb_0.4.0          parameters_0.21.7   cachem_1.1.0       
#[49] parallel_4.3.1      effectsize_0.8.8    cellranger_1.1.0    vctrs_0.6.5         jsonlite_1.8.8      hms_1.1.3          
#[55] bit64_4.0.5         Formula_1.2-5       listenv_0.9.1       systemfonts_1.1.0   jquerylib_0.1.4     parallelly_1.37.1  
#[61] glue_1.7.0          codetools_0.2-20    stringi_1.8.4       gtable_0.3.5        tables_0.9.28       munsell_0.5.1      
#[67] pillar_1.9.0        htmltools_0.5.8.1   R6_2.5.1            textshaping_0.4.0   rprojroot_2.0.4     vroom_1.6.5        
#[73] evaluate_0.24.0     highr_0.11          backports_1.5.0     memoise_2.0.1       bslib_0.7.0         Rcpp_1.0.12        
#[79] svglite_2.1.3       coda_0.19-4.1       checkmate_2.3.1     xfun_0.45           pkgconfig_2.0.3    
#
#RStudio Version: 2024.4.2.764
#R Version: R version 4.3.1 (2023-06-16)
